Ponieważ proces formowania się fali uderzeniowej w supernowej typu implozyjnego, pokazany w programi GR1D, jest bardzo skomplikowany, warto przeanalizować najprostszy znany model tego typu: równanie Burgersa:

.

Równanie to można rozwiązać np. metodą charakterystyk.

DSolve[D[u[x,t],t]+u[x,t]D[u[x,t],x]==0,u,{x,t}]

Dowolna funkcja C[1] może zostać powiązana z warunkami początkowymi, zadanymi pewną FUNKCJĄ .

DSolve[{D[u[x,t],t]+u[x,t]D[u[x,t],x]==0,u[x,0]==u0[x]},u,{x,t}]

DSolve[{D[u[x,t],t]+u[x,t]D[u[x,t],x]==0,u[x,0]==1/(1+Exp[x])},u,{x,t}]

Z braku lepszych pomysłów, jako prości fizycy bierzemy f. Gaussa.

u0=Function[x,1/(1+x^2)]

Manipulate[ContourPlot[u==u0[x - u*t],{x,-12,12},{u,0,1},AspectRatio->1/4,PlotPoints->128,ImageSize->600],{{t,0},-5,5}]

Jak widać, rozwiązanie analityczne w pewnym momencie ,,przełamuje się’’. Wykres przestaje być funkcją, i z matematycznego punktu widzenia pochodna staje się nieskończona. Formalnie dalszej ewolucji nie ma. Czas końca ewolucji można obliczyć różniczkując f. uwikłaną (Zad.).

-1/u0'[x] /.Solve[(D[u0[x],x,x]//Simplify) ==0,x]

tMAX=-1/MinValue[D[u0[x],x],x∈Reals]

Block[{t=tMAX},ContourPlot[u==u0[x - u*t],{x,-12,12},{u,0,1},AspectRatio->1/4,PlotPoints->128]]

Ale co właściwie się dzieje dla t>tMAX ? W przypadku wielu nieliniowych równań falowych (na wodzie, OTW)  pojawiają się ,,osobliwości’’, ,,wybuch’’, ,,solitony’’ itp. Aby to zrozumieć, podejmiemy naiwną próbę rozwiązania numerycznego.

rozwNUM0=NDSolve[{D[u[x,t],t]+u[x,t]D[u[x,t],x]==0,u[x,0]==u0[x],u[-12,t]==u[12,t]},u,{x,-12,12},{t,0,12tMAX}][[1]]

U0=u/.rozwNUM0

Manipulate[Plot[U0[x,t],{x,-12,12},PlotRange->{0,2},AspectRatio->1/8,ImageSize->600],{t,0,12 tMAX}]

Powyższy wynik jest w dosyć oczywisty sposób błędny. W momencie zbliżania się do ,,fali uderzeniowej’’ pojawiają się gwałtowne oscylacje numeryczne, a sam wynik jest nonsensowny. Pojawia się naturalne pytanie, czy my w ogóle rozwiązujemy prawidłowe równanie?  Z równaniem, jak się okaże, wszystko jest w porządku, ale jedna ze standardowych metod jego rozwiązywania polega na jego modyfikacji poprzez tzw. ,,sztuczną lepkość’’:

ϵ=0.002 (* Sztuczna lepkość ϵ->0 *)

Po dodaniu sztucznej lepkości równanie  zmienia charakter:

.

Nadzieja w tym, że teraz równanie rozwiążemy numerycznie bez problemu, a granica rozwiązanie dla ϵ->0 okaże się dobrze określona.

rozwNUM1=NDSolve[{D[u[x,t],t]+u[x,t]D[u[x,t],x]==ϵ D[u[x,t],x,x],u[x,0]==u0[x],u[-12,t]==u[12,t]},u,{x,-12,12},{t,0,12tMAX}][[1]]

U1=u/.rozwNUM1

Manipulate[Plot[U1[x,t],{x,-12,12},PlotRange->{0,2},AspectRatio->1/8,ImageSize->600],{t,0,12 tMAX}]

Podejmiemy teraz próbę analitycznego rozwiązanie równania Burgersa z lepkością. Co zaskakujące, jest to wykonalne.

Clear[ϵ]

burgers=D[u[t,x],t]+u[t,x]D[u[t,x],x]==0

DSolve[burgers,u,{t,x}]

burgersViscous=D[u[t,x],t]+u[t,x]D[u[t,x],x]-ϵ D[u[t,x],x,x]

travellingWave=DSolve[burgersViscous==0,u,{t,x}]

Otrzymaliśmy rozwiązanie szczególne tzw. travelling wave.

U2=u/.travellingWave[[1]]

Clear[uL,uR]

$Assumptions={ϵ>0,C[2]>0,uL>uR}

eqR=uR==Limit[U2[0,x],x->Infinity]

eqL=uL==Limit[U2[0,x],x->-Infinity]

eqL//FullForm

Reduce[{eqL,eqR},{C[1],C[2]}]

Refine[%,ϵ>0&&uL>uR]

Solve[{eqL,eqR},{C[1],C[2]},Reals]

coeffRozw=Refine[%,uL>uR]//Simplify

U2[0,0]/.coeffRozw//Simplify

U2[t,x]/.coeffRozw[[1]]/.C[3]->0//FullSimplify

U3=Function[{t,x},Evaluate[%]]

Na podstawie rozw. analitycznego możemy wyznaczyć prędkość poruszania się fali. Rozważmy dowolny punkt na fali, np. średnią (Subscript[u, L]) i (Subscript[u, R]):

U3[t,x]-(uL+uR)/2//Simplify

Ponieważ (Subscript[u, L])>(Subscript[u, R]), wyrażenie w nawiasie musi być zero. Oznacza to, że:

t (uL+uR)-2 x==0

Solve[%,x]

Jak widzimy, prędkość fali to ewidentnie .

burgersViscous/.u->U3//Simplify (* Spr. rozwiązania *)

Teraz zwizualizujemy ruch fali.

uL=2;
uR=1;

U3[t,x]

Rozwiązanie zachowuje się tak jak oczekiwalibyśmy w przypadku fali uderzeniowej.

Manipulate[Plot[(1/2) (3+Tanh[((3 t-2 x)/(8 ϵ))]),{x,-12,12},PlotRange->{0,2},AspectRatio->1/8,ImageSize->600,PerformanceGoal->"Quality",GridLines->{{(uL+uR)/2*t},None}],{t,0,12 tMAX},{{ϵ,0.01},0,0.1,0.001}]

Nawiązując do powyższego rozwiązania, możemy sformułować kluczowy w teorii  fal uderzeniowych problem Riemanna:

Jak wygląda rozwiązanie zagadnienia początkowego równania  Burgersa (w ogólności: dowolnego tego typu równania lub układu równań różniczkowych cząstkowych), gdy dane początkowe wyglądają następująco:

Manipulate[Plot[uL+(uR-uL)HeavisideTheta[x],{x,-2,2},PlotRange->2{-1,1},ExclusionsStyle->Red],{{uL,3/2},-2,2},{{uR,1/2},-2,2}]

Rozwiązanie problemu Riemanna okazuje się w ogólności bardzo nietrywialne. W przypadku r. Burgersa, dla uL>uR mamy ,,falę uderzeniową’’ poruszającą się w prawo z prędkością s=(uL+uR)/2. A co stanie się, gdy uL<uR?

burgersViscous

burgersViscous/.u->Function[{t,x},(1+Exp[(2x-t)/4/ϵ]*(1+Erf[x/Sqrt[4 ϵ t]])/(1-Erf[(x-t)/Sqrt[4 ϵ t]]))^-1]

FullSimplify[%]

r=Function[{t,x,ϵ},(1+Exp[(2x-t)/4/ϵ]*(1+Erf[x/Sqrt[4 ϵ t]])/(1-Erf[(x-t)/Sqrt[4 ϵ t]]))^-1];

Manipulate[Plot[r[t,x,ϵ],{x,-12,12},PlotRange->{0,2},AspectRatio->1/8,ImageSize->600,PerformanceGoal->"Quality"],{t,0,12 tMAX},{{ϵ,0.1},0,0.1,0.001}]

Aby problem całkowicie  wyjaśnić, portzebujemy ogólnego rozwiązania r. Burgersa. Pomimo ,,ewidentej’’ nieliniowości, okazuje się, że da się je sprowadzić do liniowego równania rozchodzenia się ciepła! Potrzebujemy to tego transformacji Cole-Hopfa.

burgersViscous/.u->Function[{t,x},-2ϵ D[Log[ϕ[t,x]],x]]//Simplify

Integrate[%,x]

Simplify[%]

Solve[%==0,Derivative[1,0][ϕ][t,x]]

heatEq=Derivative[1,0][ϕ][t,x]==(Derivative[1,0][ϕ][t,x]/.%[[1]])

DSolve[heatEq,ϕ,{t,x}]

DSolve[{heatEq,ϕ[0,x]==DiracDelta[x]},ϕ,{t,x}]

Równanie rozchodzenia się ciepła można rozw. np. za pomocą transformaty Fouriera.

rozwHeatFourier=DSolve[{heatEq,ϕ[0,x]==ϕ0[x]},ϕ,{t,x}][[1]]

Zapiszmy powyższe w nieco bardziej zjadliwej formie, z użyciem standardowych oznaczeń:

ϕ[x,t]/.rozwHeatFourier /. K[1]->ω

FourierTransform[Exp[(-ϵ ω^2+ψ[t]) t],ω,x]

Integrate[%,{x,-Infinity,Infinity}]

Plot[UnitStep[1-x]*UnitStep[x+1],{x,-12,12}]

Integrate[((E^(-((((x-y))^2)/(4 t ϵ))))/(2 (Sqrt[π]) (Sqrt[t ϵ])))*UnitStep[1-y]*UnitStep[y+1],{y,-Infinity,Infinity}]

Block[{ϵ=1,t=0.001},Plot[(((Sqrt[t]) (Sqrt[ϵ]) (Erf[((1-x)/(2 (Sqrt[t]) (Sqrt[ϵ])))]+Erf[((1+x)/(2 (Sqrt[t]) (Sqrt[ϵ])))]))/(2 (Sqrt[t ϵ]))),{x,-12,12},PlotRange->{0,1}]]

ϵ Derivative[0,2][ϕ][t,x]-Derivative[1,0][ϕ][t,x]/.ϕ->Function[{t,x},(( (Erf[((1-x)/(2 (Sqrt[t]) (Sqrt[ϵ])))]+Erf[((1+x)/(2 (Sqrt[t]) (Sqrt[ϵ])))]))/(2 ))]//Simplify

ArcTanh[3]//N

D[Log[(2/(1+(E^x)))],x]

Plot[-((E^x)/(1+(E^x))),{x,-12,12}]

Plot[x+Log[(2/(1+(E^x)))],{x,-6,6}]//N

Integrate[1/(1+Exp[x]),{x,0,Glaisher}]/.Glaisher->x

Exp[((x+Log[(2/(1+(E^x)))])/(-2ϵ))]

Block[{ϵ=1,t=1,x=0},Integrate[((E^(-((((x-y))^2)/(4 t ϵ))))/(2 (Sqrt[π]) (Sqrt[t ϵ])))*Exp[((y+Log[(2/(1+(E^y)))])/(-2ϵ))],{y,-Infinity,Infinity}]]

Block[{ϵ=1/2,t=1,x=0},Integrate[((E^(-((((x-y))^2)/(4 t ϵ))))/(2 (Sqrt[π]) (Sqrt[t ϵ])))*Exp[((y+Log[(2/(1+(E^y)))])/(-2ϵ))],{y,-Infinity,Infinity}]]

Block[{ϵ=1/2,t=1},Integrate[((E^(-((((x-y))^2)/(4 t ϵ))))/(2 (Sqrt[π]) (Sqrt[t ϵ])))*Exp[((y+Log[(2/(1+(E^y)))])/(-2ϵ))],{y,-Infinity,Infinity}]]

Block[{ϵ=1/2},Integrate[((E^(-((((x-y))^2)/(4 t ϵ))))/(2 (Sqrt[π]) (Sqrt[t ϵ])))*Exp[((y+Log[(2/(1+(E^y)))])/(-2ϵ))],{y,-Infinity,Infinity}]]

-2ϵ D[Log[(1/2) (E^(-x)) ((E^(t/2))+(E^x)) ],x]//Simplify

Manipulate[Plot[{((2 (E^(t/2)) 1/2)/((E^(t/2))+(E^x))),1/(1+Exp[x])} ,{x,-12,12}],{t,0,12,1}]

burgersViscous/.u->Function[{t,x},((2 (E^(t/2)) ϵ)/((E^(t/2))+(E^x)))]//FullSimplify

%/.ϵ->1/2

Simplify[%]

Rozwiązanie równania Burgersa metodą Godunowa, odporną na formowanie się fal uderzeniowych.

Clear[U,num]

num=64

h=1/num

xTBL=Table[(i-1/2)*h,{i,1,num}]

initPulse=Function[x,If[x≤1/4||x≥3/4,0,1]];

(*U=Table[initPulse[(i)*h],{i,1,num}]//N*)

U=Table[NIntegrate[initPulse[x],{x,(i-1)*h,i*h}]/h,{i,1,num}]//N

Total[U]*h-NIntegrate[initPulse[x],{x,0,1},WorkingPrecision->32]

ListPlot[U]

CFL=0.95

Δt=h*CFL

RotateRight[{a,b,c,d,e}]

anim1=NestList[#-Δt/h(#^2/2-RotateRight[#]^2/2)&,U,256];

anim2=NestList[#-Δt/4/h(RotateRight[#]+3*#-5*RotateLeft[#]+RotateLeft[#,2])-Δt^2/4/h^2(RotateRight[#]-#-RotateLeft[#]+RotateLeft[#,2])&,U,256];

Manipulate[
Show[
ListPlot[{xTBL,anim1[[k]]}//Transpose,PlotRange->{-1,3}],
ListPlot[{xTBL-h*0.5,anim1[[k]]}//Transpose,InterpolationOrder->0,Joined->True,PlotRange->All],
ContourPlot[u==initPulse[x-u*(k-1)*Δt],{x,0,1},{u,0,1},ContourStyle->Red,PerformanceGoal->"Quality"]]
,{k,1,Length[anim1],1}]

h*anim2[[112]]//Total

Integrate[initPulse[x],{x,0,1}]